Boltzmann and hydrodynamic description for self-propelled particles 
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We study analytically the emergence of spontaneous collective motion within large bidimensional 
groups of self-propelled particles with noisy local interactions, a schematic model for assemblies of 
biological organisms. As a central result, we derive from the individual dynamics the hydrodynamic 
equations for the density and velocity fields, thus giving a microscopic foundation to the phenomeno- 
logical equations used in previous approaches. A homogeneous spontaneous motion emerges below 
a transition line in the noise-density plane. Yet, this state is shown to be unstable against spatial 
perturbations, suggesting that more complicated structures should eventually appear. 
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Collective motion of self-propelled interacting agents 
has become in recent years an important topic of interest 
for statistical physicists. Phenomena ranging from ani- 
mal flocks (e.g. fish schools or bird flocks) to bacteria 
colonies human crowds |3|, molecular motors jij, or 
even interacting robots 0, depend only on a few gen- 
eral properties of the interacting agents 0, Q ■ From a 
physicist viewpoint, it is thus of primary importance to 
analyze generic minimal models that could capture the 
emergence of collective motion, without entering the de- 
tails of the dynamics of each particular system. In this 
spirit, Vicsek et al. proposed a simple model 0, defined 
on a continuous plane, where "animals" are represented 
schematically as point particles with a velocity vector 
of constant magnitude. Noisy interaction rules tend to 
align the velocity of any given particle with its neighbors. 
A continuous transition from a disordered state at high 
enough noise to a state where a collective motion arises 
was found numerically Q . Recent numerical simulations 
confirmed the existence of the transition, and suggested 
that the transition may be discontinuous, with strong 
finite size effects 0, other approaches, vel ocity 

vectors have been associated with classical spins lllLll2l ; 
lattice Boltzmann models have also been proposed 

However, apart from this large amount of numerical 
data, little analytical results are available. Some coarse- 
grained descriptions of the dynamics in terms of phe- 
nomenologkal h ydr odynamic equations have been pro- 
posed 0,0,0, J on the basis of symmetry and con- 
servation laws arguments. Accordingly, the coefficients 
entering these equations have no microscopic content, 
and their dependence upon external parameters is un- 
known. Renormalization group analysis |l4| and numer- 
ical studies 'iG*! confirm the presence of a nonequilibrium 
phase transition in such systems. Still, a first-principle 
analytical approach based on the dynamics of individuals 
on a continuous space is, to our knowledge, still lacking. 
Such an approach would be desirable to gain a better 
understanding of the spontaneous symmetry breaking in 
two-dimensional systems with continuous rotational sym- 



metry, a phenomenon that cannot occur in equilibrium 
systems due to the presence of long wavelength modes, 
as shown by Mermin and Wagner |rg|. Indeed, although 
the Mermin- Wagner theorem does not hold in nonequilib- 
rium system, one may wonder whether long wavelength 
modes still play an important role [l^ . 

In this short note, we introduce a microscopic bidimen- 
sional model of self-propelled particles with noisy and 
local interaction rules tending to align the velocities of 
the particles. We derive analytically hydrodynamic equa- 
tions for the density and velocity fields, within a Boltz- 
mann approach. The obtained equations are consistent 
with previous phenomenological proposals 0, 0, E, ■ 
Most importantly, we obtain explicit expressions for the 
coefficients of these equations as a function of the micro- 
scopic parameters. This allows us to analyze the phase 
diagram of the model in the noise-density plane. 

Definition of the model- We consider self-propelled 
point-like particles moving on a continuous plane, with a 
velocity vector v of fixed magnitude vq (to be chosen as 
the velocity unit) in a reference frame -hence, Galilean 
invariance no longer holds. The velocity of the particles 
is simply defined by the angle 9 between v and an ar- 
bitrary reference direction. Particles evolve ballistically 
until they experience either a self-diffusion event (a ran- 
dom "kick" ) , or a binary collision that tends to align the 
velocities of the two particles. To be more specific, the 
velocity angle 9 of any particle is changed with a prob- 
ability A per unit time to a value 9' — 9 + r/ [Fig. H^a)], 
where 77 is a Gaussian noise with distribution po(^/) smd 
variance CTq . In addition, binary collisions occur when the 
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FIG. 1: Schematic view of the dynamics of the model: (a) 
self-diffusion events, (b) binary collisions with alignment in- 
teractions -see text for notations. 
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FIG. 2: Numerical simulations with multiple-particle inter- 
actions, from IC] (dashed line in (a), (b) and (d)) and with 
binary collisions (full line), (a) average order parameter; (b) 
Binder cumulant; (c) distribution of the order parameter for 
binary collisions; (d) density profile along the coordinate x±. 



distance between two particles becomes less than do (in 
the following, we set do = The velocity angles di and 
6*2 of the two particles are then changed into 9[ = 9 + rji 
and e'2^0 + r]2 [Fig-db)], where = Arg(e*''i + e*^^) is 
the average angle, and rji and r/2 are independent Gaus- 
sian noises with the same distribution p{ri) and variance 
cr^, that may differ from ctq. 

Binary versus multiple-particle interactions - To con- 
firm that binary collisions are sufficient to capture the 
phenomena reported in numerical simulations |a, Hof , we 
performed numerical simulations of a model with binary 
collisions , and compared them with results obtained 
in a model with multi-particle interactions [l^ . In both 
models, N particles evolve on a periodic domain of lin- 
ear size L, with the same density po — N / L"^ = i in 
natural microscopic units {L = 256 for the binary col- 
lisions model, and L = 512 for the other model). The 
order parameter (f), where ip — N^^\J2j=i^^^^\: is 
shown on Fig. [2ta) as a function of the reduced noise 
e = {at — a) /at, a'^ being the variance of the noise, and 
at the value of a at the transition. Fig. El^b) shows the 
Binder cumulant G = 1 — ((/5*)/3((^^)^. The negative peak 
indicates a discontinuous transition toward spontaneous 
motion, which is confirmed in Fig. HJc) by plotting the 
probability distribution function (PDF) of the order pa- 
rameter (for binary collisions) for e below, above and very 
close to the transition [i^l ■ The distribution is clearly bi- 
modal at the transition, which is typical of discontinuous 
transitions. Finally, Fig.[2Id) presents the density profile 



obtained when spontaneous motion sets in, indicating the 
presence of a stripe with higher density. This stripe is es- 
sentially invariant along the direction, and is moving 
along the direction (on which the profile is measured, 
using a moving frame). Note that the profile is asymmet- 
ric, with a higher slope on the front. Thus, a model with 
only binary collisions is legitimate and behaves qualita- 
tively in a similar way as a model with more complicated 
interactions. 

Boltzmann equation- The Boltzmann equation de- 
scribing the evolution of the one-particle phase-space dis- 
tribution f{r,9,t) reads 



dl 
dt 



(r, 0, t) +e{e)- V/(r, 9, t) - I^lf] + /eoi [/] (1) 



where /dif[/] accounts for the self-diffusion phenomenon, 
and /coi[/] describes the effect of collisions; e{9) is the 
unit vector in the direction 9. /dif[/] is given by 



/dif[/] = -A/(r,0,t) + A / dO' dr^poiv) (2) 
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The collision term /coi[/] is evaluated as follows. By def- 
inition, two particles collide if their relative distance be- 
comes less than do- In the referential of particle 1, par- 
ticle 2 has a velocity V2 = e(92) — 6(6*1). Thus, particles 
that collide with particle 1 between t and t-\-dt are those 
that lie, at time t, in a rectangle of length jvjl and of 
width 2do. This leads to 



IcoAf] ^ -f{r,9,t) d9'\e{9')-e{9)\f{r,9',t) (3) 
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with 9 = Arg(e**i -f e**^). It can be checked easily that 
the miiform distribution f(r,9,t) = p/^n, is a solution 
of Eq. 1^ for any density, and whatever the form of the 
noise distributions po(j]) ^'^'^ P{v)- 

Hydrodynamic equations - Let us now define the hy- 
drodynamic density and velocity fields p(r, t) and u(r, t) 



p{v,t) = / d9f{r,9,t) (4) 

J —TT 

p(r,t)u(r,i) = / d9f{v,9,t)e{9) (5) 



Integrating the Boltzmann equation over 9 yields the 
continuity equation for p(v, t) 



dp 
'di 



V • (pu) = 



(6) 
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The derivation of a hydrodynamic equation for the ve- 
locity field is less straightforward, and involves an ap- 
proximation scheme. Let us introduce the Fourier series 
expansion of /(r, 0, t) with respect to 9 



fkir,t) 



d0f{r,e,t) e 



ike 



(7) 



Multiplying Eq. ^ by e^''^ and integrating over 9 leads to 
an infinite set of coupled equations for /fc(r,i). We note 
that, identifying complex numbers with two-dimensional 
vectors so that e*^ corresponds to e{9), the Fourier co- 
efficient fi{r,t) is nothing but the "momentum" field 
w(r,t) = p{r,t) u{r,t). Thus the evolution equation 
for /i(r,t) should yield the hydrodynamic equation for 
u{r,t). Yet, as fki^.t) is coupled to all others fi{r,t), a 
closure relation has to be found. In the following, we 
assume that the velocity distribution f{r,9,t) is only 
slightly non-isotropic, or in other words that u(r, t) is 
small as compared to the individual velocity of particles, 
and that the hydrodynamic fields vary on length scales 
that are much larger than the microscopic length do. As 
a result, the velocity equation is obtained from the equa- 
tion for /i through an expansion to leading orders in 
/fe(r, i) and in space and time derivatives. Noting that 
/o(r,t) = p{r,t) = 0{1), we set /i(r,i) = 0(e), e « 1. A 
consistent scaling ansatz, confirmed by a numerical inte- 
gration of Eq. Q in the steady state, is /fc(r, t) = C'(el'^l). 
Expanding to order e^, one only keeps the terms in /i and 
/2 in the evolution equation for fi. A similar expansion 
of the equation for /2 leads to a closure relation for the 
equation on /i, finally leading to the following hydrody- 
namic equation ^| 



dw 
'dt 



■ 7(w • V)w 



1 



(8) 



+ {l^i — ^w^)w + !/V^w — k(V • w)w 



where the different coefficients are given by 
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The first term in the r.h.s. of Eq. may be thought of 
as a pressure gradient, introducing an effective pressure 
p = \{p — KvP). The second term describes the local 
relaxation of w, whereas the third term corresponds to 




FIG. 3: Phase diagram of the model in the plane (p/A, a). A 
transition line (full line: ao = cr; dashed line: ao — 1) indi- 
cates the linear instability threshold of the state w — |w| — 0. 



the usual viscous term, and the last one may be inter- 
preted as a feedback from the compressibility of the flow. 
The fact that 7^1 (apart from special values of cr) 
in the advection term expresses that the problem is not 
Galilean invariant. Note that v, 7 and k are always pos- 
itive; /i can change sign, and ^ > whenever fi > 0. 
All the terms are compatible with the phenomenological 
equation of motion of Toner et al. |l4| . However, our 
approach provides explicit forms for the coefficients. In 
particular, the coefficient in front of the term V(V • w) 
is strictly zero in our case. Besides, there is no term of 
the form (w • V)^w due to the order of truncation of the 
Boltzmann equation. 

Phase diagram.- We can now study the spontaneous 
onset of collective motion in the present model. As a 
first step, it is interesting to consider possible instabil- 
ities of the spatially homogeneous system, that is the 
appearance of a uniform nonzero field w. Equating all 
space derivatives to zero leads to the simple equation 



'dt 



(/i — ^w^)w 



(14) 



Clearly, w = is solution for all values of the coefficients, 
but it becomes unstable for /i > 0, when a nonzero solu- 
tion Wq — y/jjj^ e appears, where e is an arbitrary unit 
vector. From Eq. ifT^ . the value p = corresponds to a 
threshold value pt 



Pt 



7rA(l 



-'^0/21 



4(e^ 



V2 



(15) 



The transition line defined by pt in the plane (p/A, a) 
is plotted on Fig. O for (Tq = a and for a fixed value 
do = 1. If (Tq — a, the instability occurs at any density, 
provided the noise is low enough. On the contrary, at 
fixed ao, the instability disappears below a finite density 
P( = 37rA(l — e~'^o/^)/4. Both transition lines saturate 
at a value at = (2 In §)i/2 w 0.90. 
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Let us now test the stability against perturbations of 
the above spatially homogeneous flow w(r, t) = wq, with 
finite density po, in an infinite space. From Eq. (|14|l . 
it is clear that Wq is stable against spatially homoge- 
neous perturbations. Yet, this solution may be unstable 
against finite wavelength perturbations \22\. To check 
this issue, we introduce a perturbation around the ho- 
mogeneous steady-state solution 

p{r,t) ^ po + 5p{r,t), w(r, i) = Wo + (5w(r,t) (16) 

and linearize Eq. ^ in 6p{r,t) and dw{r,t). Linear sta- 
bility is then tested with the ansatz 

6p{r, t) = 6po e^*+*'^■^ (5w(r, t) = Swo e^*+*'^■^ (17) 

where q is a given wavevector, by looking for the dis- 
persion relation s(q). By choosing 6wq and q along the 
same direction as wg, one finds for the real part of s 

for small |q|, with hq = 4(e~'^^/^ — |)/7r and wq = \wq\, 
indicating the onset of a long wavelength instability since 
3?(s) becomes positive at small enough |q| The 
spatially homogeneous states w = and w = wq are 
thus both unstable, so that more complicated structures 
should eventually appear in the system. The "stripes" of 
higher density moving over a low density background, re- 
ported in [lOj, may be examples of such patterns. Phys- 
ically, the instability may be interpreted as follows: if 
locally p > pq [p < pq), the local velocity u increases (de- 
creases), creating velocity gradients that generate a den- 
sity wave. Note that the perturbations that destabilize 
the long-range order correspond to longitudinal waves, 
at odds with what happens in the two-dimensional XY- 
model |l9l | which might be thought of as an equilibrium 
counterpart of the present model ■ 

Discussion.- Our analytical approach has several ad- 
vantages when compared with pure numerical simula- 
tions of similar microscopic models. First, the hydro- 
dynamic equations may be used to get analytical solu- 
tions in reference cases with simple geometries, and to 
analyze their stability against perturbations. Second, in 
more complicated situations, these equations may be in- 
tegrated numerically, allowing one to study much larger 
systems than with direct simulations of the particles. 

The hydrodynamic equations (|6I8(I have been derived 
from a Boltzmann approach and their validity is in prin- 
ciple restricted to a low density regime (note that the 
transition may occur at low density by choosing A <C 1). 
However, as verified for many systems, the validity of the 
Boltzmann equation often goes well-beyond the a priori 
expected limit. One also expects that in this low density 
regime, the hydrodynamic equations should not depend 
strongly on the details of the interactions, as the shape of 



p{ri). Another limitation comes from the assumption that 
w is small. This assumption is valid to describe the evolu- 
tion of small perturbations around the zero velocity state. 
When crossing the transition line, the assumption is self- 
consistent, as the resulting homogeneous velocity field 
grows continuously from zero. Yet, the finite- wavelength 
instability may drive the system into a regime where the 
approximation is no longer valid. Checking this point 
requires to find the structures emerging from Eqs. (|6I8|I . 
and to compare them with numerical simulations. Work 
in this direction is under investigation (l^. 
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